Artificial quantum-dot Helium molecules: 
Electronic spectra, spin structures, and Heisenberg clusters 
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Energy spectra and spin configurations of a system of N = 4 electrons in lateral double quantum 
dots (quantum dot Helium molecules) are investigated using exact diagonalization (EXD), as a 
function of interdot separation, applied magnetic field (B), and strength of interelectron repulsion. 
As a function of the magnetic field, the energy spectra exhibit a low-energy band consisting of a 
group of six states, with the number six being a consequence of the conservation of the total spin 
and the ensuing spin degeneracies for four electrons. The energies of the six states appear to cross 
at a single value of the magnetic field, and with increasing Coulomb repulsion they tend to become 
degenerate, with a well defined energy gap separating them from the higher-in-energy excited states. 
The appearance of the low-energy band is a consequence of the formation of a Wigner supermolecule, 
with the four electrons (two in each dot) being localized at the vertices of a rectangle. Using spin- 
resolved pair-correlation distributions, a method for mapping the complicated EXD many-body wave 
functions onto simpler spin functions associated with a system of four localized spins is introduced. 
Detailed interpretation of the EXD spin functions and EXD spectra associated with the low-energy 
band via a 4-site Heisenberg cluster (with B-dependent exchange integrals) is demonstrated. Aspects 
of spin entanglement, referring to the well known A-qubit Dicke states, are also discussed. 

PACS numbers: 73.21. La, 31.15.V-, 03.67.Mn, 03.65.Ud 



I. INTRODUCTION 

The field of two-dimensional (2D) semiconductor 
quantum dots (QDs) has witnessed rapid expansion 
in the last several years, both experimentallyis 2 - and 
theoretically^^ 5 - Along with fundamental interest in the 
properties of such systems, and as a test ground for highly 
correlated electrons, a major motivation for these grow- 
ing endeavors has been the promising outlook and po- 
tential of quantum dots concerning the implementation 
of solid-state quantum computing and quantum informa- 
tion devices ,SiL§i2ii£ To this effect highly precise control 
of the space and spin degrees of freedom of a small num- 
ber N of confined electrons (down to an empt y 11 ' 12 ' 13 
QD) needs to be achieved, and experimentally this was 
demonstrated recently for two electrons in a lateral dou- 
ble quantum dot molecule (see Ref. [H, and references 
therein). From the theoretical standpoint, high-level 
computational methods that reach beyond the level of 
mean-field approximation are needed; 5 - with the ability 
to provide solutions that preserve all the symmetries of 
the many-body Hamiltonian, and in particular those as- 
sociated with the total spin. In this context, electrons in 
quantum dots exhibit localization in space and formation 
of Wigner molecules (see, e.g., Ref.[H). When the spin de- 
gree of freedom is considered, such Wigner molecules may 
be viewed as finite Heisenberg spin cluster a 14 ' 15 whose 
quantum behavior (due to finite-size fluctuations and 
correlation effects) differs drastically from the behavior 
expected from magnetic systems in the thermodynamic 
limit 

There is an abundance of experimental and theoreti- 
cal publications concerning circular single quantum dots 
with a small number of electrons ,LiLiaS a 17 i 18 i 19 In this 



paper, we use exact diagonalizatio n 4 ' 5 ' 18 (EXD) to in- 
vestigate the properties of lateral double quantum dots 
(DQDs) containing four electrons. DQDs are referred to 
also as artificial molecules. Specifically in the case of four 
electrons they can be viewed as artificial quantum dot He- 
lium molecules £Q DQDs containing two electrons have 
been already studied extensively both experimentally 2 
and theoretically . 5 ' 21 ' 22 However, experimental studies 
of DQDs with more than two electrons are relatively 
fewj 23 ' 24 We are aware of a single theoretical EXD study 
of a lateral DQD with three electrons^ and another 
one of two laterally coupled quantum rings with three 
electrons^ 

In light of the novel quantum behavior discovered in 
our investigations (compared to circular QDs with re- 
gard to the spectra, spin structures, and analogies with 
Heisenberg clusters), we hope that the present work 
would serve as an impetus for further experimental stud- 
ies on lateral DQDs. In particular, as a function of the 
magnetic field, we find that: (1) A low-energy band of six 
states develops as the strength of the Coulomb repulsion 
increases, separated by an energy gap from the other ex- 
cited states, and (2) All six states appear to "cross" at 
a single value of the magnetic field. The crossing point 
gets sharper for larger interdot distances. We find that 
the specific number of crossing states (six) derives from 
the spin degeneracies and multiplicities expressed by the 
branching diagram i2I 

The formation of the low-energy band is a consequence 
of the localization of the four electrons within each dot 
(with two electrons in each dot). This localization leads 
to formation (with increasing strength of the Coulomb 
repulsion) of a Wigner supermolecule^ 2 ^ with the four lo- 
calized electrons (two in each dot) being located at the 
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corners of a rectangular parallelogram (RP). Using the 
spin-resolved pair-correlation functions, we show how to 
map the EXD many-body wave functions onto the spin 
functions associated with the four localized spins. This 
mapping leads us naturally to study analogies with finite 
systems described by a model Heisenberg Hamiltonian 
(often referred to as finite Heisenberg clusters). Specif- 
ically, we provide a detailed interpretation of the EXD 
spin functions and EXD spectra associated with the low- 
energy band via a 4-site finite Heisenberg cluster charac- 
terized by two (intradot and interdot) exchange integrals. 

More importantly, our EXD calculations exhibit a 
prominent oscillatory magnetic-field dependence of the 
two exchange integrals entering in this 4-site Heisen- 
berg Hamiltonian. Such strong independence of the ex- 
change integrals has been found in previous theoretical 
studies in the simpler case of two electrons in double 
quantum dots^ ' 21 i 22 ' 29 i 30 as well as in anisotropic sin- 
gle quantum dots i 31 ' 32 The independence in the case of 
two electrons in quantum dots has also been observed 
exp er iment ally 3 1 ' 3 2 1 33 Following an earlier proposal, 7 the 
B-dependence of the exchange integral for the two- 
electron case has developed into a central theme in ex- 
perimental efforts aiming at solid-state implementation 
of quantum computing. Our EXD results in this pa- 
per extend the independence of the exchange integrals 
to larger numbers (N > 2) of electrons in quantum dot 
molecules. 

We further discuss that the determination of the equiv- 
alent spin functions enables consideration of aspects of 
entanglement regarding the EXD solutions. In par- 
ticular, we show that the formation of Wigner super- 
molecules leads to strongly entangled states known in 
the literature of quantum information as 7V-qubit Dicke 
states^i 3 ^! 3 ! 

We finally mention that the trends in the excitation 
spectra (e.g., formation of a low-energy band) and en- 
tanglement characteristics (e.g., mapping to spin func- 
tions of localized electrons) found in the case of double 
quantum dots have many analogies with those found in 
other deformed configurations, and in particular single 
anisotropic quantum dots; see, e.g., the case of three elec- 
trons in Ref. [H. 

The plan of the paper is as follows: 

• Section II describes the two-dimensional two- 
center-oscillator (TCO) external confining poten- 
tial that models the double quantum dot. 

• Section III reviews the many-body Hamiltonian 
and the exact-diagonalization method as imple- 
mented in this paper. 

• Section IV outlines some theoretical background 
regarding the general form of four-electron spin 
eigenfunctions and the branching diagram (which 
describes the break-down of spin multiplicities for 
given N). 
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FIG. 1: (Color online) Cuts at y = of the two-center- 
oscillator confining potential along the x-axis for an inter- 
dot separation of d — 60 nm. Solid line (blue): double- 
oscillator confining potential with a smooth connecting neck 
specified by e b = 14/ Vo = 0.5. Dashed line (black): double- 
oscillator confining potential without a smooth connecting 
neck. Dotted line (red): double-oscillator confining poten- 
tial with a smooth connecting neck specified by e b — 6. The 
remaining parameters are hu) x i = Tiu) X 2 = huio = 5.1 meV, 
m* = 0.070m e , hx=hi = 0, Vbi = V02 = V , e\=e\ = e h . 

• Section V describes our numerical results from the 
exact diagonalization, that is, the EXD spectra, the 
electron densities, and the spin-resolved conditional 
probability distributions (CPDs). 

• Section VI provides an interpretation of the numer- 
ical EXD results for the 6-state lower-energy band 
with the help of a 4-site Heisenberg Hamiltonian. 

• Section VII contains two parts which discuss (a) 
the importance of _B-dependent exchange integrals 
in the Heisenberg Hamiltonian and (b) the aspects 
of entanglement exhibited by the EXD wave func- 
tions. 

• Section VIII offers a summary. 

• Finally, the Appendix describes the single-particle 
spectrum of the two-center oscillator (when the 
two-body Coulomb interaction is omitted; non- 
interacting model). 



II. TWO-DIMENSIONAL 
TWO-CENTER-OSCILLATOR CONFINING 
POTENTIAL 

In the two-dimensional two-center-oscillator, the 
single-particle levels associated with the confining po- 
tential of the artificial molecule are determined by the 
single-particle Hamiltonian 39 

tt m , 1*22,1*2/2 

H = T + -m uj y y + -m uj xk x k 

+ V neck (x) + h k + ^-B ■ s , (1) 
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where x' k — x — Xk with k = 1 for x < (left) and 
k = 2 for x > (right), and the /i^'s control the relative 
well-depth, thus allowing studies of hetero-QDMs. y de- 
notes the coordinate perpendicular to the interdot axis 
(x). T = (p - eA/c) 2 /2m*, with A = 0.5(-By, Bx, 0), 
and the last term in Eq. (fTJ) is the Zccman interaction 
with g* being the effective g factor, /ib the Bohr magne- 
ton, and s the spin of an individual electron. The most 
general shapes described by H are two semiellipses con- 
nected by a smooth neck V nec k(x) (see solid line in Fig. 
[T]). xi < and x 2 > are the centers of these semiel- 
lipses, d = x 2 — x\ is the interdot distance, and to* is the 
effective electron mass. 

For the smooth connecting neck, we use V nec k{x) = 
\m*ul k [C k x'% + V k x' k 4 }6(\x\ - \x k \), where 9{u) = for 
u > and 9(u) — 1 for u < 0. The four constants C k 
and T> k can be expressed via two parameters, as follows: 
C fe = (2 - Ae b k )/x k and V k = (1 - 3e b k )/x 2 k , where the 
barrier-control parameters e\ = (Vb — h k )/Vok are related 
to the actual (controlable) height of the bare interdot 
barrier (V&) between the two QDs, and Vok — TO * w xfc :E fe/2 
(for hi = h 2 , V i = V 02 = Vo). 

The single-particle levels of H, including an external 
perpendicular magnetic field B, are obtained by numer- 
ical diagonalization in a (variable-with-separation) basis 
consisting of the eigenstates of the auxiliary (zero-field) 
Hamiltonian: 

n 2 1 1 

= ^7 + 2 u yV + 2 m UxkXk + k ■ ( ' 

The eigenvalue problem associated with the auxiliary 
Hamiltonian [Eq. @] is separable in x and y, i.e., the 
wave functions arc written as 



<Pi(x,y) = X^(x)Y n (y), 



(3) 



with i = {/i, n}, i = 1, 2, . . . , if . 

The y n (y) are the eigenfunctions of a one-dimensional 
oscillator, and the X^(x < 0) or X f _ l (x > 0) can be 
expressed through the parabolic cylinder function a 40 ! 41 
Ui'fk, (-l) fc ^ fc ], where £ fc = x' k y/2m*uj xk /h, j k = (-E x + 
h k ) / (tiijj xk ) , and E x = (p, + 0. b)%uj x i + hi denotes the x- 
eigenvalues. The matching conditions at x = for the 
left and right domains yield the x-eigenvalues and the 
eigenfunctions X fl (x). The n indices are integer. The 
number of fi indices is finite; however, they are in gen- 
eral real numbers. 

In the Appendix, we discuss briefly the energy spec- 
tra associated with the single-particle states of the two- 
center oscillator Hamiltonian given by Eq. (TT|) . We follow 
there the notation presented first in Ref. |42|. For further 
details, see Ref. 

In this paper, we will limit ourselves to QDMs with 
x 2 = —Xi and tkuy — hu x i — huj x2 = hujQ. How- 
ever, in several instances we will compare with the case 
of a single elliptic QD where x 2 — —x\ — and 
hujy =/= hcu x = huj x i = huj x2 . In all cases, we will use 
huj = 5.1 meV, to* = 0.070m e (this effective-mass value 
corresponds to GaAs), and K = 50 (which guarantees 
numerical convergence* 4 -) . 



III. THE MANY-BODY HAMILTONIAN AND 
THE EXACT DIAGONALIZATION METHOD 

The many-body Hamiltonian 7i for a dimeric QDM 
comprising N electrons can be expressed as a sum of the 
single-particle part H(i) defined in Eq. |T]) and the two- 
particle intcrclcctron Coulomb repulsion, 



N 



N N 9 



(4) 



i=l j>i 



where k is the dielectric constant and denotes the 
relative distance between the i and j electrons. 

As we mentioned in the introduction, we will use the 
method of exact diagonalization for determining 4 ^ the so- 
lution of the many-body problem specified by the Hamil- 
tonian ([I]). 

In the EXD method, one writes the many-body wave 
function $^ XD (ri, r 2 , . . . , r N ) as a linear superposition 
of Slater determinants \& (rj., r 2 , . . . , rjy) that span the 
many-body Hilbert space and are constructed out of the 
single-particle spin-orbitals 



Xj [x,y) = <pj (x, y)a, if 1 < j < K , 



and 



Xj(x,y) = (pj- K (x,y)P, if K < j < 2K, 
where a(/3) denote up (down) spins. Namely 



(5) 



(6) 



4> 



EXD 

N,q 



(ri 



r N ) = y]cj^(n,...,r N ), (7) 



where 



N 



Xji( r i) •■■ Xj«(>i) 



Xh(r N ) ... Xj N ( r N) 



(8) 



and the master index I counts the number of arrange- 
ments {ji, j 2 , . . • , Jn} under the restriction that 1 < ji < 
32 < ■ ■ ■ < 3n < 2K. Of course, q = 1,2,... counts the 
excitation spectrum, with q — 1 corresponding to the 
ground state. 

The exact diagonalization of the many-body 
Schrodinger equation 



(9) 



transforms into a matrix diagonalizatiom problem, which 
yields the coefficients Cj and the EXD eigenenergies 
i?]v* D ■ Because the resulting matrix is sparse, we im- 
plement its numerical diagonalization employing the well 
known ARPACK solver^ 

The matrix elements (^ r ] v -|W|^ , ^ r ) between the basis 
determinants [see Eq. ©] are calculated using the Slater 
rules4£ Naturally, an important ingredient in this respect 
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are the two-body matrix elements of the Coulomb inter- 
action, 



dr 1 dr 2 ip*{r 1 )^*(r 2 )- 



1*1 - r 2 



(10) 

in the basis formed out of the single-particle spatial 
orbitals <Pi(r), i = 1,2,..., K [Eq. (0)]. In our ap- 
proach, these matrix elements are determined numeri- 
cally and stored separately. Varying the dielectric con- 
stant k and/or the interdot-barrier parameter e b does not 
require a recalculation of the Coulomb- interaction matrix 
elements, as long as the remaining parameters are kept 
the same. 

The Slater determinants [see Eq. ([SJ] conserve the 
third projection S z , but not the square S 2 of the total 
spin. However, because S 2 commutes with the many- 
body Hamiltonian, the EXD solutions are automatically 
eigenstates of S 2 with eigenvalues S(S +1). After the 
diagonalization, these eigenvalues are determined by ap- 

=,EXD 

-N,q 



plying S 2 onto * D and using the relation 



S 2 *f 



(N a -N p ) 2 /4 + N/2 + Y, 



*<3 



(ii) 



where the operator Wij interchanges the spins of electrons 
i and j provided that their spins are different; N a and Np 
denote the number of spin- up and spin-down electrons, 
respectively. 

Of great help in reducing the size of the matrices to 
be diagonalized is the fact that the parity (with respect 
to the origin) of the EXD many-body wave function is 
a good quantum number for all values of the magnetic 
field when h\ — h 2 - Specifically, the rcy-parity operator 
associated with reflections about the origin of the axes is 
defined as 



x iJ ^ X g D (ri,r2,r3,r 4 ) 



V xll <5> 



(-ri,-r 2 ,-r 3 ,- 



t 4 ) 
(12) 

and has eigenvalues ±1. 

One can also consider partial parity operators V x and 
V v associated solely with reflections about the x and y 
axis, respectively; of course V xy = V X V V . We note that 
unlike V xy , the partial parities V x and V v are conserved 
only for zero magnetic fields (B = 0). With the two- 
center oscillator cartesian basis that we use [see Eq. ([5)1], 
it is easy to calculate the parity eigenvalues for the Slater 
determinants, Eq. ([8]), that span the many-body Hilbcrt 
space. Because X^ix) and Y n (y) conserve the partial V x 
and Vy parities, respectively, one finds: 



V X y^ = (-)£*=! "H+»«$W 



(13) 



where m, and n% count the number of single-particle 
states associated with the bare two-center oscillator [see 
the auxiliary Hamiltonian Hq in Eq. ([2"])] along the x axis 
and the simple oscillator along the y direction (with the 




FIG. 2: The branching diagram for the spin degeneracies. 
The total-spin quantum number S is given on the vertical 
axis, and the number of particles, N, on the horizontal one. 
The numbers inside the circles give the number, g(N,S), of 
linear independent (and orthogonal) spin functions for the 
corresponding values of N and S. 



assumption that the lowest states have m — and n = 0, 
since they are even states). We note again that the in- 
dex \i in Eq. ([3]) is not an integer in general, while m 
here is indeed an integer (since it counts the number of 
single-particle states along the x direction). 



IV. MANY-BODY SPIN EIGENFUNCTIONS 

For completeness and for the reader's convenience, we 
outline in this section several well established (but of- 
ten not well known) properties of the many-body spin 
eigenfunctions which are useful for analyzing the trends 
and behavior of the spin multiplicities exhibited by the 
EXD wave functions for N = 4 electrons. We stress 
here that the ability to describe spin multiplicities is an 
advantage of the EXD method compared to the more 
familiar spin-density functional approaches whose single- 
determinantal wave functions preserve only the third pro- 
jection S z of the total spin, and thus are subject to "spin 
contamination" errors. As we will discuss below, the spin 
multiplicities of the EXD wave functions lead naturally 
to formation of highly entangled Dicke states , 34 i 35 i 36 i 37 
and most importantly to analogies with finite Heisenberg 
clusters i 14 ' 15 

A basic property of spin eigenfunctions is that they ex- 
hibit degeneracies for N > 2, i.e., there may be more than 
one linearly independent (and orthogonal) spin functions 
that are simultaneous eigenstates of both S 2 and S z . 
These degeneracies are usually visualized by means of the 
branching diagram 27 displayed in Fig. [2j The axes in this 
plot describe the number N of fermions (horizontal axis) 
and the quantum number S of the total spin (vertical 
axis). At each point (N, S), a circle is drawn containing 
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the number g(N, 5) which gives the degeneracy of spin 
states. It is found^i that 



N 



( N/2 - S 



N 

N/2 -5-1 



(14) 



Specifically for N — 4 particles, there is one spin eigen- 
function with 5 = 2, three with 5=1, and two with 
5 = 0. In general the spin part of the EXD wave func- 
tions involves a linear superposition over all the degener- 
ate spin eigenfunctions for a given 5. 

For a small number of particles, one can find compact 
expressions that encompass all possible superpositions. 
For example, for TV = 4 and 5 = 0, S z = one has4£ 



Finally, for N — 4 and 5 = 2, S z — (maximum 
polarization) case, one has: 

X20 = 

mm + lint) + mm + mm + inn) + inn) 



V6 



(17) 



V. EXACT-DIAGONALIZATION RESULTS 



A. Energy spectra 



(HI 



1 sin^| ttU> + (icos0- 



sin0)|nn) 



-(-COS0 + 



-(-cose + 



+(-cos( 



sin 0)1 nit) 

sine) mn) 
sin^ inn) - 



sin 6» mm, 

(15) 



where the parameter 8 satisfies — ir/2 < 8 < ir/2 and is 
chosen such that 8 = corresponds to the spin function 
with intermediate two-electron spin S12 = and three- 
electron spin 5i23 = 1/2; whereas 8 = ±tt/2 corresponds 
to the one with intermediate spins 5i2 = 1 and 5i23 = 
1/2. 

For N = 4 and 5 = 1, 5 Z = one has: 
X10 = 

(\j\ s[n0s[n( p - sin cos ( p~ ^ c ° s ^)inn) 

+(^5 sin6< sin - sine cos y> + i cose) |T4T4) 

+(v 77; sine cos <p - \/~ sinesintp - ^cose) ||tlt) 

V 12 V o 2 

+(v 77; sin e cos <p - \/- sine sin <p + ^cose) |tllT) 

V 12 V o 2 

+ ( \J\ sin 8 sin tp + ^ sin 6> cos tp) 1 | U 4 ) 

-(y^ sin e sin 93 + sine cos tp) |||tt) , (16) 

where the parameters 8 and y satisfy — ir/2 < 8 < tt/2 
and — 7r/2 < ip < Tt/2. Three independent spin functions 
with definite intermediate two-electron, S12, and three- 
electron, 5i23, spin values correspond to the 8 and tp 
values as follows: for 5i2 = and 5i23 = 1/2, 8 = 0; for 
5i2 = 1 and 5i 23 = 1/2, = ±tt/2 and tp = 0; and for 
5i2 = 1 and 5i23 = 3/2, 8 = ±tt/2 and tp = ±tt/2. 



The excitation spectra as a function of the applied 
magnetic field for four electrons in a double QD with in- 
terdot distance d = 2x2 = —2xi = 30 nm and no voltage 
bias between the dots [hi = I12 = 0, see Eq. ([1}] are plot- 
ted for three different values of the interelectron repulsion 
strength, i.e., weak [k = 12.5 (GaAs); see Fig. [3(a)], in- 
termediate [n = 6; see Fig. [3(b)], and strong [n = 2; see 
Fig. [3(c)] Coulomb repulsion. The interdot barrier pa- 
rameter was taken as e h = 0.5 (because hi = /12 = 0, one 
has e\ = e h 2 = e b ; see Section [U for the definitions). In 
all cases, we calculated the eight lowest energy levels. 

We observe that the lowest six levels form a band that 
separates from the rest of the spectrum through the open- 
ing of a gap. This happens already at a relatively weak 
interelectron repulsion, and it is well developed for the 
intermediate case (n = 6). It is of interest to note that 
the number of levels in the band (six) coincides with the 
total number of spin eigenfunctions for TV = 4 fermions, 
as can be seen from the branching diagram displaying 
the spin degeneracies. In particular, there is one level 
with total spin 5 = 2 (and parity V xy = 1), three levels 
with total spin 5=1 (two with V X y = 1 and one with 
Vxy = —1), and two levels with total spin 5 = (one 
with V X y = 1 and the second with V xy = — 1)- All these 
six levels approximately "cross" at one point^ situated 
at about B w 3.5 T for n = 12.5, B « 2.2 T for n = 6, 
and Brj1.1T for k = 2. 

The trends associated with the opening of a gap and 
the formation of a six-state low band appear further re- 
inforced for the larger interdot distance of d = 60 nm 
(displayed in Fig. 2] for the three values of the dielectric 
constant k = 12.5, 6, and 2, respectively). It is remark- 
able that the six lower curves "cross" now at a sharply 
defined point^ (situated at B sa 3.3 T for k = 12.5, 
B w 2.1 T for k = 6, and B « 1.0 T for k = 2). The 
six curves demonstrate additional near degeneracies re- 
grouping approximately to three curves before and after 
the crossing point, which results in a remarkable simpli- 
fication of the spectrum. 

For strong repulsion (k = 2), all six states in the low 
band are almost degenerate for both distances [d = 30 
nm; see Fig. [3(c) and d = 60 nm; see Fig. 5(c)]. 
This is a consequence of the formation of a near-rigid 
Wigner molecule (WM) with strongly localized electrons. 
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FIG. 3: (Color online) Energy spectra (as a function of the 
magnetic field B) for TV = 4 electrons in a double quantum 
dot with interdot separation d = 30 nm. (a) n — 12.5 corre- 
sponds to GaAs. (b) k = 6. (c) K = 2. The Zeeman term was 
neglected, and thus all states with the same total spin S dif- 
ferent spin projections S z are degenerate. Remaining param- 
eters: e b — 0.5, TvjJq = 5.1 meV, m* = 0.07m e , hi = hi = 0. 
For all figures in this paper, hioo, m", and hi = hi are 
the same. Energies are referenced to Nh^/ui^ + wf/4, where 
uj c — eB/(m*c) is the cyclotron frequency. 
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FIG. 4: (Color online) Energy spectra (as a function of the 
magnetic field 73) for TV = 4 electrons in a double quantum 
dot with interdot separation d — 60 nm. (a) k = 12.5 (weak 
interelectron repulsion) corresponds to GaAs. (b) k — 6. 
(c) k = 2. The interdot barrier corresponds to e b = 0.5. 
The inset in (c) displays all the lowest 12 P xy — 1 and the 
lowest 10 P xy = — 1 energies, and it demonstrates a trend 
toward formation of higher bands. Energies are referenced 
to Nhy^ + cJc/4, where uj c — eB/(m*c) is the cyclotron 
frequency. 



Namely, the overlaps between the orbitals of localized 
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electrons are very small (see, e.g., Ref . [28h . yielding small 
exchange contributions in the total energies, 2 - and thus 
all six possible spin multiplicities tend to become degen- 
erate in energy. Furthermore the physical picture of a 
near-rigid Wigner molecule suggests that the energy gap 
to the next group of states corresponds to excitation of 
the lowest stretching vibrational mode of the 4-electron 
molecule. 



Since the main panels in Figs. [3] and [3] display only the 
four lowest-in-energy states with positive parity and the 
four corresponding states with negative parity, one needs 
a larger part of the spectrum to ascertain whether higher 
bands are formed. To this end, the inset in Fig. rjjc) 
displays the twelve lowest-in-energy curves with P xy = 1 
and the ten corresponding curves with P xy = — 1. In 
addition to the 6-state low band, the inset indicates for- 
mation of a higher band comprising a total of 12 states 
(not labeled); however, a detailed study of this higher 
band falls outside the scope of the present paper. 



It is natural to anticipate at this point that the above 
behavior of the low-energy EXD spectra at low B can be 
generalized to an arbitrary number of electrons N in a 
double QD. Namely, as the strength of the interelectron 
interaction increases, a low-energy band comprising all 
possible spin multiplicities will form and it will become 
progressively well separated by an energy gap from the 
higher excitations. For example, for N = 6, an inspec- 
tion of the branching diagram in Fig. [5] leads us to the 
prediction that there will be 20 states in this low-energy 
band. A similar behavior emerges also in the case of a 
single, but strongly anisotropic quantum dot; indeed a 
low-energy band of three states (see the branching dia- 
gram in Fig. [2|) has been found for N = 3 electrons in 
Ref. M 



It is of interest to contrast the above behavior of the 
excitation spectra in a double QD with that of an N- 
electron circular dot. Specifically, in the circular QD, 
large inetelectron repulsion leads to formation of a near- 
rigid rotating Wigner molecule that exhibits a rigid mo- 
ment of inertia. Then the states inside the low-energy 
band (two states for N = 2, three for N = 3, six for 
N = 4, etc.) do not become degenerate in energy, but 
form an yrast rotational band 49 specified by L 2 /2Jo, 
where L is the total angular momentum and j7o is the 
classical moment of inertia. We note that the energy 
splittings among the yrast rotational states are much 
smaller than the vibrational energy gap in circular dots 
associated with the quantum of energy of the 

stretching (often referred to as breathing) mode of the 
polygonal-ring configuration of the quasiclassical Wigner 
molecule^^i^ 2 - 



B. Electron densities 

The electron density is the expectation value of the 
one-body operator 

N 

p(r) = ^5(r-r i ), (18) 

that is: 

p(r) = (^\p(r)\^) 

= £ CfCj(V? \p(v)M). (19) 
i, J 

Since p(r) is a one-body operator, it connects only 
Slater determinants and "J^ that differ at most by 
one spin orbital \j ( r ) j f° r the corresponding Slater rules 
for calculating matrix elements between determinants for 
one-body operators in terms of spin orbitals, see Table 
2.3 in Ref. El 

In Figs. 02a-f), we display (for the aforementioned 
three strengths of interelectron repulsion) the ground- 
state electron densities for for N — 4 electrons in the 
case of a double dot at zero magnetic field with inter- 
dot separations d = 30 nm (left column) and d — 60 nm 
(right column). 

For the weak interaction case (k — 12.5) at B — 0, the 
electron densities do not exhibit clear signatures of for- 
mation of a Wigner molecule for either interdot distance, 
d = 30 nm [Fig. Ufa)] or d = 60 nm [Fig. Mb)]. The 
Wigner molecule is well formed, however, in the case of 
the intermediate Coulomb repulsion [n = 6; see Figs.[5jc- 
d)]. One observes indeed four humps that correspond to 
the four localized electrons; they are located at (±34.88 
nm, ±13.13 nm) in the c? = 60 nm case. In the case of 
strong Coulomb repulsion (k = 2) and for the same inter- 
dot distance d = 60 nm, the electrons are further local- 
ized as can be seen from Fig.JSff); the four humps occur 
now at (±39.86 nm, ±21.02 nm). The Wigner molecule 
is also well formed in the the strong-repulsion and d — 30 
nm case, as can be seen from Fig.JSJe), with the localized 
electrons located at (±29.28 nm, ±21.11 nm). 



C. Spin-resolved conditional probability 
distributions 

1. Definitions 

In the regime corresponding to a well-defined Wigner 
molecule, the electron densities (see Sect. IV B[) are char- 
acterized by four humps that reflect the localization of 
the four electrons in the double quantum dot. Such 
charge densities do not provide any information concern- 
ing the spin structure of each EXD state. In fact, all 
six EXD states in the lower band exhibit very similar 
four-humped electron densities. 



8 




s=o 



(C) ^ K=6 (d) 



MB- 




d=30 nm 



d=60 nm 



FIG. 5: (Color online) Electron densities at B = for the 
ground state (with S = 0, S z — 0, parity P xy = 1) for iV = 4 
electrons in a double quantum dot with interdot separations 
d = 30 nm (left column) and d — 60 nm (right column) . The 
top (a-b), middle (c-d), and bottom (e-f) rows correspond 
to k = 12.5 (weak repulsion), k = 6 (intermediate repul- 
sion), and k = 2 (strong repulsion), respectively. Ground- 
state energies: (a) E = 27.609 meV, (b) E = 21.572 meV, (c) 
E = 49.217 meV, (d) E = 39.799 meV, (e) E = 111.361 meV, 
(f) E = 94.516 meV (compare Figs. EH3J. The interdot bar- 
rier corresponds to e 6 = 0.5. Distances in nm. Vertical axis 
in arbitrary units (with the same scale for all six panels). 



The spin configurations associated with a given (S, S z ) 
EXD state in the WM regime can be explored with the 
help of the spin-resolved two-point anisotropic correla- 
tion function defined as: 



P aao (r,r Q ) = 
(^ X 9 D | E S ( F ~ r ^( r ° - 'i) Wo** 1*^), (20) 

with the EXD many-body wave function given by equa- 
tion ©. 

Using a normalization constant 

Af(a,a ,r ) = J P CTCTo (r, r )dr, (21) 

we further define a related conditional probability distri- 




d=60 nm 

ground 



8 b =0.5 



d=60 nm 

ground 



8 b =6 




d=60 nm 8 =0.5 
excited 



d=30 nm 8 =0.5 
ground 



FIG. 6: (Color online) CPDs V n at B = for several EXD 
states with S = 0, S z =0, and parity P xy = 1 of N = 4 
electrons in a double quantum dot with interdot separations 
d = 60 nm (a-c) and d = 30 nm (d) . Case of strong Coulomb 
repulsion (k = 2) with an interdot barrier e b = 0.5 (a,c-d) and 
e b = 6 (b). Panels (a-b,d) correspond to ground states. Panel 
(c) corresponds to the excited second S = state for the same 
parameters as in panel (a) [see Fig. |4jc) and the branching 
diagram in Fig. [2]. Energies: (a) E — 94.516 meV, (b) E = 
96.811 meV, (c) E = 95.017 meV, and (d) E = 111.361 meV 
[compare Figs. EJc) and H£c)]. Distances in nm. Vertical 
axis in arbitrary units (with the same scale for all panels in 
Figs. [6] — [8j). The fixed point is located at the maximum of the 
hump in the lower-left quadrant of the corresponding electron 
density, i.e., at ro =(—40 nm, —21 nm) for panels (a-c) and 
ro =( — 29 nm, —19 nm) for panel (d). 



bution (CPD) as 

Vaa {r,r ) = F CTCTo (r,r )/7V(cr,cro,ro), (22) 

having the property J P CTCT0 (r, ro)dr = 1. The spin- 
resolved CPD gives the spatial probability distribution 
of finding a second electron with spin projection a under 
the condition that another electron is located (fixed) at 
ro with spin projection <7o; o and can be either up (f) 
or down Q). 

To calculate P aao (r, ro) in Eq. (j2"D|) , we use a sym- 
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metrized operator 

f CTCT0 (r,r ) = 

[<5(r-r. i )(5(ro-r J )J (T(Ti (5 (To(Tj .- 



S(r - r,)5(r - Ti)S aaj S aoa ^\ , (23) 



yielding 



P CTff0 (r,r ) = (^ x 9 D |f|^ x g D ) 

i.j 



(24) 



Since T aao (r, tq) is a two-body operator, it connects 
only Slater determinants 9f and that differ at most 
by two spin orbitals Xji ( r ) and Xh ( r )> f° r the correspond- 
ing Slater rules for calculating matrix elements between 
determinants for two-body operators in terms of spin or- 
bitals, see Table 2.4 in Ref. IM 



2. Examples of S = 0, S z = EXD states at B = 

For each charge density corresponding to a given state 
of the system, one can plot four different spin-resolved 
CPDs, i.e., V\^, V^i, V^, and Vn- This can potentially 
lead to a very large number of time consuming compu- 
tations and an excessive number of plots. For studying 
the spin structure of the S — 0, S z = states at B = 0, 
however, we found that knowledge of a single CPD, taken 
here to be Vfi (see Fig. [5]), is sufficient in the regime of 
Wigner-molecule formation. Indeed, the specific angle 
9 specifying the spin function Xqo Eq. (|15p correspond- 
ing to the CPDs portrayed in Fig. E] can be determined 
through the procedure described in the following: 

We designate with roman indices I, II, III, and IV 
the four quadrants of the (x, y) plane, starting with the 
upper left quadrant and going clockwise [see Fig. Eta)]. 
In the case of a 4e Wigner-molecule, a single electron is 
localized within each quadrant. The same roman indices 
designate also the positions of the localized electrons in 
each of the six Slater determinants (e.g., | HID, I tltl), 
etc.) that enter into the spin function A^o in Eq. (fT5"|) . We 
take always the fixed point to correspond to the fourth 
(IV) quadrant [bottom left in Fig. Et a )]- An inspection 
of Eq. (|15p shows that only three Slater determinants 
in X 00 contribute to V\i, namely | HU), I T4T4), and 
| ITTI); these are the only determinants in Eq. (|15p with 
a down spin in the 4th quadrant. From these three Slater 
determinants, only the first and the second contribute 
to the conditional probability II-f^(7) of finding another 
electron with spin-up in quadrant I; this corresponds to 
the volume under the hump of the EXD CPD in quadrant 
I [see, e.g., the hump in Fig. [5(a)]. Taking the squares 
of the coefficients of | Hil) and I Titi) in Eq. (fT5|) . one 
gets 



1 



• cos 6 



sin 9 



(25) 



Similarly, one finds that only | TT44) and | ITT4) con- 
tribute to H^i(II), and that 



n n (jj) cx 



sm 



I 



■ cost 



sin( 



(26) 



Integrating under the humps of the EXD CPD in 
quadrants I and //, we determine numerically the ratio 
IL[l(I)/Hfl(II), which allows us to specify the absolute 
value of 9 (within the interval -90° < 9 < 90°) via the 
expressions in Eqs. (|25|) and (f26|) . The restriction to the 
absolute value of 9 is a result of the squares of the sine 
and cosine entering in II-pj_ (I) and (II). To obtain the 
actual sign of 9, additional information is needed: for ex- 
ample the ratio Hn(I)/Il^i(III) can be used in a similar 
way, where 



n n (/J7) cx ( I cost 



1 

- cost 
2 



-sin0 



sm t 



(27) 



Using the method described above, we find that 9 i=s 
—60° for the EXD ground state at d = 60 nm (larger 
interdot distance) and n = 2 [strong repulsion; sec Fig. 
[SJa)], and the corresponding spin function simplifies to 



4o = -§l TtU> + §l UU) + §I 4UT)-§| UTT>. (28) 

Remarkably, increasing the interdot barrier from e b = 
0.5 [Fig. Eta)] to e b = 6 [Fig. Efb)], while keeping the 
other parameters constant, does not influence much the 
composition of the associated spin function, which re- 
mains that given by Eq. (|28p . This happens in spite 
of the visible change in the degree of localization in the 
electronic orbitals, with the higher interdot-barrier case 
exhibiting a sharper localization. 

In Fig. Etc), we display the Vn CPD for an excited 
state with S = 0, S z = (having P xy = I and energy 
E = 95.017 meV), with the remaining parameters being 
the same as in Fig. Et a )- For this case, following an 
analysis as described above, we found the angle 9 « 30°, 
which is associated with a spin function of the form 



1 



1 



1 



1 



y{2) _ 



^imD + ^ium-^itut) 



1 



2V3 



UH) 



(29) 



We note that the spin functions in Eqs. (|28|) and (|29|) 
are orthogonal. 

In Fig. Etd), we display the V-\\ CPD for the ground 
state with S = 0, S z =0 (having P xy = 1 and energy 
E = 111.361 meV) and for the shorter interdot distance 
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d=60 nm 

(b) ! 




and energy E = 94.757 meV, at the larger interdot 
separation d = 60 nm. Again we consider the case of 
strong Coulomb repulsion (k = 2) with an interdot bar- 
rier e b = 0.5. The corresponding spin function Xiq [Eq. 
(fT6|) ] depends on two different angles 9 and cf>, and one 
needs at least two different CPDs for determining their 
specific values. For this purpose, we display also the Vn 
CPD for the same state in Fig.[7fb). 

The specific values of 9 and 0" associated with the CPDs 
in Figs. [7[a) and [Tfb) can be determined through the 
ratios n n (7)/n u (77) and II n (7)/II n (777) [associated 
with Fig. Ha)] andn u (/)/n u (/7) and n u (7)/n u (777) 
[associated with Fig.[7Jh)], where 



n u (7) oc ism 2 6»sin 2 
3 



12 



■ sin 2 9 cos 2 c 



cos 2 9 



V2 



■ sm f sin cos < 



V6 



sm t) cos f sm ( 



12 



sm v cos a cos < 



(31) 



FIG. 7: (Color online) CPDs at B = for excited EXD states 
with S = 1, S z = 0, and parity P xy = 1 of N = 4 electrons 
in a double quantum dot at the larger interdot separation 
d = 60 nm (a-b) and the shorter interdot separation d — 30 
nm (c). Panels (a) and (c) display V\i CPDS (down-up), 
while panel (b) displays a different Vn CPD (down-down), 
but for the same state as in (a). Case of strong Coulomb 
repulsion (k = 2) with interdot barrier e b — 0.5. Energies: (a- 
b)E = 94.757 meV, and (c) E = 111.438 meV [compare Figs. 
[3jc) and[4jc)]. Distances in nm. Vertical axis in arbitrary 
units (with the same scale for all panels in Figs. [6] — |5J. The 
fixed point is located at the maximum of the hump in the 
lower- left quadrant of the corresponding electron density, i.e., 
at ro =(—40 nm, —21 nm) for panels (a-b) and ro =( — 29 nm, 
— 19 nm) for panel (c). Note that this is a case with S = 1; 
the previous Fig. [6] displayed S = cases. 



d = 30 nm. For this case, we found an angle 9 w —63.08° 
which corresponds to the following spin function: 



X$ = -0.5148| TTU>+ 0.4838| UU> + 0.031| TUT) 
0.031| ITU) + 0.4838| HIT) - 0.5148| ||TT). 

(30) 

From a comparison of the above result with that for 
the larger d = 60 nm [see Eq. ([28]) ]. we conclude that the 
difference in interdot distance results in a slight variation 
of the spin functions. 



3. Examples of S = 1, S z = EXD states at B = 

In this section, we turn our attention to partially po- 
larized EXD states with 5 = 1. 

In Fig. [TJa), we display the V n CPD at B = for 
an excited state with S = 1,S Z = 0, parity P xy = 1, 



15 1 

Pit I (II) oc - sin 2 9 sin 2 6 H sin 2 9 cos 2 d>+ - cos 2 i 



v/2 . 



1 



■ sm f sm cos < 



H — == sm # cos 9 cos <z>, 
V12 



V6 



sm f cos W sm ( 



(32) 



U n (III) ex - sin 2 sin 2 + - sin 2 9 cos 2 . 
3 6 



\pi i 

■— — sin 2 9 sin 1 cos H — cos 2 (33) 



and 

n u (/) cx 



(34) 



n il( // ) « | \/ - sm « mum - 



sm y cos < 



■ cos 9 



(35) 



U U (III) cx I sin # sin + sin cos j . (36) 



Using Eqs. ((31J) — (|36|) and the numerical val- 
ues of the ratios H n (I)/U n (H) and n n (/)/n Ti (m) 
n u (7)/n u (77) and n u (7)/n u (777) (specified via a 
volume integration under the humps of the EXD CPDs), 
we determined that 9 — —45° and sin^ = — -\/2/3, 
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S=0 B=2T 




lower energy higher energy 



FIG. 8: (Color online) V n CPDs for the EXD states at B = 2 
T with S = 0, S z = 0, and parity P xy = 1 of N = 4 electrons 
in a double quantum dot at the larger interdot separation 
d = 60 nm. (a) The lower energy of the two S — states 
(see branching diagram in Fig. [2} . (b) Higher energy S — 
state. Case of strong Coulomb repulsion (n = 2) with interdot 
barrier e b = 0.5. Energies: (a) E = 94.605 meV and (b) E = 
95.047 meV [compare Fig. 0(c)] . Distances in nm. Vertical 
axis in arbitrary units (with the same scale for all panels in 
Figs. [6] — [8]). The fixed point is located at ro =(—40 nm, —21 
nm) . 



cos0 = vO/3 (i.e., <f> w -54.736°). Thus, the corre- 
sponding spin function reduces to the simple form 

*io = UU> - y|l iUT). (37) 

In Fig. [7(c), we display the P u CPD at S = for a 
similar excited state as in Fig. [7(a) (with S — 1, 5 2 = 0, 
parity P xy = 1, and energy E = 111.438 meV) of N = 4 
electrons at the shorter interdot separation d = 30 nm. 
Here too we consider the case of strong Coulomb repul- 
sion (k = 2) with interdot barrier e b = 0.5. We note that 
the localization of electrons is stronger for the larger in- 
terdot distance [compare Fig. [7(a) with Fig. [7(c)]. This 
difference, however, does not influence the coefficients en- 
tering into the associated spin function, which we found 
to remain very close to the specific form in Eq. (|37[) . 



4- Spin-resolved conditional probability distributions at 
B=£Q 

In Fig. [5] we display EXD CPDs at a finite value of 
the magnetic field, and precisely at B = 2 T, for the 
two states of the low-energy band with S = 0, S x = 
(at the larger interdot separation d = 60 nm and strong 
interelectron repulsion n = 2). This value of B was cho- 
sen to lie beyond the crossing point for the six states 
of the low-energy band [which happens at B ~ 1 T; see 
Fig. 0(c)] . Comparison with the CPDs of the correspond- 
ing states at zero magnetic field [see Figs.[B(a) and[H(c)] 
shows that the spin structure of the associated Wigner 



molecule varies rather slowly with the increasing mag- 
netic field in the range < B < 2.5 T. 

Following the height of the humps in the left upper 
quadrants, one observes that the CPD in Fig. [5(a) (case 
of lower- energy state at B — 2 T with S = and P xy = 1) 
corresponds to that of Fig. [H(a) (case of lower- energy 
state at B = with S = and P xy = 1). Similarly, 
the CPD in Fig. [5(b) at B = 2 T {higher- energy state) 
corresponds to that of Fig. [6(c) at B — (higher- energy 
state). From these results, we concludes that the two 
states with S = and P xy = 1 do not really cross at 
the 'crossing' point at B ~ 1 T. In reality, this point is 
an anticrossing point for these two states, although the 
anticrossing gap is too small to be seen with the naked 
eye. This behavior agrees with that expected from states 
having the same quantum numbers. We checked that a 
similar observation applies for the two other states in the 
low-energy band having the same quantum numbers, i.e., 
those having 5 = 1 and P xy = — 1. 

VI. ANALOGIES WITH A 4-SITE 
HEISENBERG SPIN CLUSTER 

In Section IV CI using the spin-resolved CPDs, we 
showed that the EXD many-body wave functions in the 
Wigner-molecule regime can be expressed as a linear su- 
perposition of a small number of Slater determinants and 
that this superposition exhibits the structure expected 
from the theory of many-body spin functions. This find- 
ing naturally suggests a strong analogy with the field of 
nanomagnets and quantum magnetism, usually studied 
via the explicitly spin-dependent model effective Hamil- 
tonian known as the Heisenberg Hamiltonian ) 14 ' 15 i 53 
given by: 7i' H = Ylij JijSi'Sj, where are the ex- 
change integrals between spins on sites i and j. Even in 
its more familiar, simpler form 

Hh = ^ JySj-Sj, (38) 

(hj) 

that is that of the spin- 1/2 Heisenberg antiferromagnet 
with nearest-neighbor interactions only and Jy > 0, it is 
well known that the zero-temperature (at B = 0) solu- 
tions of Hamiltonian (|38p involve radically different forms 
as a function of the geometry, dimensionality, and size. 

Generalizing this behavior to finite magnetic fields B, 
we have found that the rich variety of the EXD energy 
spectra presented in Figs. [5] and [H as well the EXD spin 
functions of Section IV CI can be related to those of a 4- 
site Heisenberg Hamiltonian H^f (B) with B-dependent 
exchange constants Jij(B) = JijFij(B), and with the 
four electrons being located at the vertices of a rectan- 
gular parallelogram (RP) as discussed earlier. Due to 
the reflection symmetry, Tl^f(B) has only two different 
exchange constants J\2 = J34 and J14 = J23, i.e., 

nl F (B) = j 12 (s)(s r s 2 + s 3 -s 4 ) + 

Ju(£)(SrS 4 + S 2 -S 3 ), (39) 
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where 1 -> /, 2 II, 3 //J, and 4 -> TV [in a clock- 
wise direction, see Fig. HJa)]. Since the B = spin ex- 
change interaction constants are expected to decrease 
exponentially 54 with the distance between the two sites i 
and j, one expects that the Heisenberg model will repro- 
duce the present EXD results in the regime Jn « Ju- 
To proceed, it is sufficient to use the six-dimensional 



Ising Hilbert subspace for zero total-spin projection 
(S z = 0), which is spanned by the following set of 
basis states [we follow here the ordering in Eq. (fT"5|) ]: 

|i) - 1 tni), 12) - 1 tin), |3) - 1 tut), K) - 1 im>, 

|5> I ITIT), and 1 6) -> | ||tt>- In this subspace, the 
Heisenberg Hamiltonian given by Eq. (|39[) ] can be writ- 
ten in matrix form as 



H*?{B) = 





/ 


J\2 — Jli 


Jli 








Ju 





\ 






Jl4 


-{JX2 + Jli) 









Jl4 




1 







Jl2 


Jl4 — J12 





Jl2 







2 







Jl2 







Jl2 











Jli 





Jl2 




-(Jl2_+ Jw) 


Jl4 






V 





Jli 








Jl4 


•^12 — </l4 


/ 



r 



(40) 



A lengthy, but straightforward, calculation yields the 
general eigenvalues £ i and corresponding eigenvectors Vj 
of the matrix (|40| . The eigenvalues are: 



£i = -(Ji4 + Ji 2 )/2, 



(41) 



v 6 = {i,-y,-i + y,-i + y,-y,i}, s = o, (53) 

where 

* = r + Q(l,r), (54) 



£2 = {Ju - Ji 2 )/2, 



£3 = (J12 - Ju)/2, 



£4 = (Jl4 + Jl2)/2, 



£e = -(Jw + .A 2 )/2 + Q(Ji4, Ju), 



where 



Q{a, b) = \Ja 2 -ab + b 2 . 



(42) 



(43) 



(44) 



y = r-Q(l,r), 



(55) 



£s = -{Ju + Ji2)/2-Q{JuJ l2 ), (45) 



(46) 



(47) 



The corresponding unnormalized eigenvectors and 
their total spins are given by: 

Vi = {0,-1,0,0,1,0}, 5 = 1, (48) 



V 2 = {0, 0,-1, 1,0,0}, 5=1, (49) 



V 3 = {-1,0,0,0,0,1}, 5 = 1, (50) 



V 4 = {1,1, 1,1, 1,1}, 5 = 2, (51) 



V 5 = {l,-X,-1 + X,-1 + X,-X,l}, 5 = 0, (52) 



and r = Jyij Ju- 

To understand how the Heisenberg Hamiltonian in Eq. 
(f3"9"| captures the rich behavior seen in the EXD spectra 
of Figs. [3] and [4] we start with the limiting case r — > 0, 
which is applicable (see below) to the larger interdot dis- 
tance d = 60 nm. In this limit, one can neglect J12 
compared with Ju , which results in partial degeneracies 
within the band; namely one has £2 = £i = £q — Ju/2, 
£1 = £ 3 = - J u /2, and £ 5 = -3 J u /2. This 3 - 2 - 1 de- 
generacy pattern is independent of the magnetic-field de- 
pendence through Fi4{B) [jjy — JijFij{B)] and is char- 
acteristic of all three EXD spectra (for k = 12.5, 6, and 
2) associated with the larger interdot distance d = 60 
nm. Furthermore, the fact that all six curves in the 
EXD lowest-energy band appear to cross at the same 
point B c (k) (reversing at the same time the order of the 
degenerate levels) suggests that 



F lA {B) ~ cos[7rB/(2Sf 4 )]. 



(56) 



It is remarkable that the behavior described above is 
prominent even for the weak interelectron repulsion n = 
12.5 [see Fig. HJa)] when the extent of electron local- 
ization and the formation of a Wigner molecule are not 
clearly visible via an inspection of the corresponding elec- 
tron densities [see Fig. [5Jb)]. 

Of interest is the fact that the ability of the Heisenberg 
Hamiltonian in Eq. (f3"5)) to reproduce the EXD trends is 
not restricted solely to energy spectra, but extends to the 
EXD wave functions as well. Indeed when J12 — > 0, the 
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last two eigenvectors of the Heisenberg matrix (having 
5 = 0) become 

V 5 -►{1,-1,0,0, -1,1}, (57) 

and 

V 6 -{l,l,-2,-2,l,l}. (58) 

When multiplied by the normalization factor, the wave 
functions represented by the eigenvectors in Eqs. (|57|) and 
([S"5]) coincide (within an overall =f1 sign) with the EXD 

spin functions X^ and <Y y in Eqs. (|2"5]1 and ([2"9"]l . re- 
spectively. In addition, when again multiplied by the cor- 
responding normalization factor, the wave function rep- 
resented by the eigenvector Vi [Eq. ijlgjl] (having total 
spin 5 = 1 ) coincides (within an overall minus sign) with 
the EXD spin function X 10 in Eq. ([37]) . 

The EXD spectra and spin functions for the shorter 
distance d — 30 nm can be analyzed within the frame- 
work of the 4-site Heisenberg Hamiltonian (|40p when 
small (compared with J14), but nonnegligible, values of 
the second exchange integral J12 are considered. In this 
case, the partial three-fold and two-fold degeneracies are 
lifted. Indeed in Figs. H a) (k = 12.5) and^b) (k = 6), 
the EXD lowest-energy band consists of six distinct lev- 
els. For the strong interelectron case with k = 2 and 
d = 30 nm [Fig. EJc)], however, the EXD spectra indi- 
cate that the effective (selfconsistent-field) potential bar- 
rier between the dots due to the Coulomb repulsion is 
high enough to reduce J12 to a negligeable value and to 
produce spectra exhibiting the characteristic 3 — 2 — 1 
degeneracy pattern that is prominent in the spectra as- 
sociated with the larger interdot distance d = 60 nm. 
(For the wave functions in the k — 2 and d = 30 nm 
case, however, the influence of J12 cannot be neglected, 
see below.) 

In addition to the lifting of the partial degeneracies 
within the lowest-energy band, one observes from an ex- 
amination of Figs. [3£a) and [31(b) the occurrence of two 
characteristic secondary oscillations (as a function of B) 
emerging out of the previously three-fold and two-fold de- 
generate levels. Using the Heisenberg Hamiltonian, these 
secondary oscillations can be described by taking the sec- 
ond exchange constant to have a i?-dependence similar 
to that in Eq. (|56"j). i.e., 

F 12 (B) ~ caa[*B/{2B&) + 0o], (59) 

with < Bi 4 , 4>q being a phase shift. This secondary 
oscillation is superimposed on the main oscillation spec- 
ified by F\±(B) [Eq. ([56]) ] in accordance with the expres- 
sions for the Heisenberg energy levels £2 [Eq. (|4"2"| ]. £4 
[Eq. gl], £ 6 [Eq. (USD], and £ Y [Eq. @TJ], £ 3 [Eq. O]. 

Another characteristic feature developing for d = 30 
nm in Figs. [3fa) and^b) is the anticrossing gap A be- 
tween the two 5 = states. According to the Heisenberg 
model this gap is given by A = 2 J 12 cos[7ri?^/(2i3^) + 
4>q] . As a concrete example of the above, we estimated 



that the spectrum of Fig. [3jb) can be rather well re- 
produced using the expressions for the Heisenberg eigen- 
values when J12/J14 ~ -1/4.1, B^ 2 sa B^/2.5, and 

00 = tt/2.4. 

Similar to the findings in the larger-distance (d = 60 
nm) case, the agreement between EXD and Heisenberg- 
model results in the smaller d = 30 nm distance includes 
also the wave functions. Indeed, as discussed in Section 
IVC31 the P xy = 1, S = 1 EXD spin function for d = 30 
nm (and B = 0, n = 2) was found to be identical to the 
one determined for the larger interdot distance d = 60 
nm, i.e., it is given by expression (|3"7| . This reflects the 
remarkable property that the Heisenberg eigenvector Vi 
[Eq. (|48p ] is independent of the two exchange constants 
J14 and J12, and thus independent of the interdot dis- 
tance d (as well as of the dielectric constant k and the 
magnetic field B). On the other hand, the 5 — Heisen- 
berg eigenvectors [Eqs. ([52]) and ([53]) ] do depend on the 
ratio r = J12/J14, which is in agreement with the fact 
that the EXD ground-state spin function Xqq for d = 30 
nm in Eq. (|30|) is slightly different from the correspond- 
ing 5 — EXD spin function X^ 1 for d = 60 nm [see 
Eq. ([28]) ]. With consideration of the normalization fac- 
tor, we estimate that the Heisenberg eigenvector V5 [Eq. 

([S"2")) ] agrees with X^ 1 when r ~ —1/7.5. 

It is of interest to contrast the EXD spin functions 
determined in Section IV CI with the well known solu- 
tions of the Heisenberg Hamiltonian [Eq. ([39])] when 
the four spins are located on four sites arranged in a 
perfect square ; 15 ' 16 i.e., when the two exchange inte- 
grals are equal; J14 = J12 = J > 0. (The perfect- 
square arrangement arises^ 5 - also in the case of forma- 
tion of a four-electron Wigner molecule in a single cir- 
cular quantum dot.) In this case, the ground state of 
Hh is the celebrated resonating valence bond (RVB) 
stat o 15 ' 16 which forms the basic unit block in many theo- 
retical approaches aiming at describing high-temperature 
superconductors^ 6 - The RVB state has quantum numbers 
5 = 0, 5 2 =0 and a Heisenberg energy £ 5 = —2 J; it is 
given by the normalized version of V5 [Eq. (|52p ] when 
r = 1 , that is by (see also Refs. [TH and [la ) 

^ VB = 27I lmi> ~7I ITm> + 2v! ITUT> 

(60) 

Although the (excited-state) EXD X$ [Eq. ([29])] in 
the quantum-double-dot case portrayed in Fig. [B](c) ap- 
pears (superficially) to be similar to the (ground-state) 
RVB A^ VB [Eq. f5D]l]. the two are not equal. Indeed the 
coefficients of the pair of Slater dctcrmimants | tit J.) an d 

1 Itlt) have been interchanged with those of | |ttl) an d 
I tilt)- Similar observations apply also to the remaining 
pair of 5 = and S z = states that are orthogonal to 
Xq^ [see Xqq in Eq. ([2"5]k case of double quantum dot] 
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and to A' ( ^ VB (case of a perfect square), respectively. 

We note that the differences in the Abo spin functions 
between the DQD case (corresponding to a rectangular 
parallelogram) and the perfect-square case are also re- 
flected in the V n CPDs. Indeed the CPDs of the DQD 
(Fig. [5]) exhibit equal-height humps along the smaller 
side of the parallelogram, while those of the perfect- 
square configuration (and/or circular quantum dot) ex- 
hibit equal-height humps along a diagonal. 55 

VII. DISCUSSION 

A. Magnetic-field dependence and relevance to 
quantum computing 

Strongly correlated electrons on a lattice are frequently 
described by the Hubbard-model Hamiltonian 

^Hubbard = _ £ t . . c | ctC . ct + V J- ni]nil , (61) 

where is the hopping integral from site j to site i, 
a = ±1 [equivalently this denotes a spin up (|) or a 
spin down (J,)], and n ia — c\ a c ic! , with cJ CT and c i(J being 
single-particle creation and annihilation operators for the 
site i. U is the on-site Coulomb repulsion. 

It is well known that the one-band Hubbard model at 
half-filling reduce a 57 i 58 (to lowest order) to a Heisenberg 
antifcrromagnetic Hamiltonian Hh [see Eq. (1551) ] in the 
limit of the on-site Coulomb repulsion U > being large 
relative to the hopping integral Uj . In the absence of an 
applied magnetic field (B — 0), t%j can be taken to be 
real and the corresponding exchange integrals are given 
byS 

jHubbard = 2tl/U, (62) 

In the presence of a magnetic field (which is the 
case of this paper), ty picks up 59 i 60 a Peierls phase 
exp[(ze/c?i) A A • dr], where A is the vector potential. 
In this case is complex, and one must replac e 59 ! 60 
tfj — > hjt*j in Eq. (f6"2"| . The complex conjugation, how- 
ever, cancels any magnetic-field effect associated with the 
Peierls phase factor, which means that the jH ubbard are 
independent of B. 

In sharp contrast with this Hubbard-model result, our 
EXD calculations indicate that the exchange integrals 
entering in the Heisenberg Hamiltonian of Eq. (f55| de- 
pend strongly on the magnetic field. Such strong in- 
dependence of the exchange integrals has been found in 
previous theoretical studies in the simpler case of two 
electrons in double quantum dots^ ' 21 i 22 ' 29 i 30 as well as 
in anisotropic single quantum dotS i 31 i 32 (For two elec- 
trons, the exchange integral is calculated as the energy 
difference between the singlet and triplet states.) This in- 
dependence in the case of two electrons in quantum dots 
has also been observed experimentally. 31 ! 32 ! 33 Following 



an earlier proposal, 7 the independence of J for the two- 
electron case has developed into a central theme in exper- 
imental efforts focussing on solid-state implementation of 
quantum computing j2i24 In this context, our EXD results 
in this paper extend the -B-dependence of the exchange 
integrals to larger numbers (TV > 2) of electrons in quan- 
tum dot molecules. 

The physics underlying the emergence of such strong 
B-dependence in the case of solid-state artificial nanos- 
tructures is clearly related to the importance^ of orbital 
magnetic effects resulting from the much larger size (by 
a factor of 10000) of the electronic wave functions in 2D 
quantum dots compared to that in natural atoms. For 
spin interactions between electrons localized within the 
natural atoms, huge magnetic fields (of order 10000 T) 
are required for reproducing a B-dcpcndcncc of the ex- 
change integrals similar to that discussed in this paper. 



B. Aspects of spin entanglement 

Neel antiferromagnetic ordering, where the average 
spin per site < S? >= ( — l) J+1 /2, is an important mag- 
netic phenomenon in the thermodymanic limit^S associ- 
ated with breaking of the total-spin symmetry. The finite 
size magnetic clusters discussed here exhibit a sharply 
different behavior in this respect. Indeed, as discussed 
in Ref. [l6j, the four-site Neel state is the single Slater 
determinant | 4T4T) ( or T4T4))- It is clear that the total- 
spin conserving EXD functions Abo (as well as the corre- 
sponding Heisenberg eigenvectors) are multideterminen- 
tal and have an average spin per localized electron (per 
site) < S z j >= 0. 

We concur with Ref. [l6| that the phenomenon of Neel 
antiferomagnetism is radically modified in assemblies of 
few electrons. In this section, we argue that instead 
of "antiferromagnetic ordering" the appropriate physical 
concept for the WM states found earlier is that of spin en- 
tanglement. Indeed, in the previous sections, we showed 
that the EXD wave functions in the regime of Wigner- 
molecule formation can be approximated as a superpo- 
sition of a small number of Slater determinants corre- 
sponding to well structured spin functions; see, e.g., Xqq 
in Eq. (|2"5)) . This is a great simplification compared to 
the initial EXD superposition [Eq. J?}], where the count- 
ing index is usually I > 500, 000. This reduction of the 
molecular EXD solutions to their equivalent spin func- 
tions (described in Section IV CP (or to the Heisenberg 
eigenvectors described in Section IVI[) enables one to in- 
vestigate their properties regarding fundamental quan- 
tum behavior associated with quantum correlations and 
fluctuations beyond the mean field. 

The mathematical theory of entanglement is still devel- 
oping and includes several directions. One way to study 
entanglement is through the use of properly defined mea- 
sures of entanglement, e.g., the von Neumann entropy 
which utilizes the single-particle density matrix. Another 
way is to catalog and specify classes of entangled states 
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that share common properties regarding multipartite en- 
tanglement. A well known class of A-qubit entangled 
states are the Dicke states , 34 ' 35 i 36 ' 37 which most often 
are taken to have the symmetric form: 

X%t={1) (| 11^1 00... 0) + Perm). (63) 

Each qubit is a linear superposition of two single-particle 
states denoted by or 1, and the symbol 'Perm' stands 
for all remaining permutations. The or 1 do not have 
to be necessarily up or down 1/2-spin states. Two-level 
atoms in linear ultracold traps have already been used 
as an implementation of a qubit. Dicke states appear 
in many physical processes like superradiance and super- 
fluorescence. They can also be realized with photons, 
where the qubits correspond to the polarization degree 
of freedom^ 

In the 1/2-spin case of fermions (e.g., for electrons), the 
Dicke states of Eq. (|63|) correspond to a fully symmetric 
flip of k out of N localized spins. It is apparent that the 
four-qubit fully polarized (S = 2 with spin projection 
S z =0) EXD solution is reproduced by X 2 o of Eq. ([17]). 
and thus it is of the symmetric Dicke form (with k = 2) 
displayed above in Eq. (|6"3"|) . On the other part, the DQD 
EXD states (with S z = 0) studied in Section IV CI with 
S = and/or 5 = 1 represent a natural generalization of 
Eq. (|63p to the class of asymmetric Dicke states. 

Dicke states with a single flip (k = 1) arc known as 
W stateSi 61 i 62 For N = 4 electrons, the latter states are 
related to EXD solutions with S z = ±1. For the con- 
nection between W states and EXD states for N = 3 
electrons in anisotropic quantum dots, see Ref. l38l . W 
states have already been realized experimentally using 
two-level ultracold ions in linear traps^ 3 - 

VIII. SUMMARY 

Extensive investigations of lateral double quantum 
dots containing four electrons (artificial quantum-dot 
Helium molecules) were performed using the exact- 
diagonalization method (described in Section |TTT| , as a 
function of interdot separation, applied magnetic field, 
and strength of intcrclcctron repulsion. Novel quantum 
behavior was discovered compared to circular QDs con- 
cerning energy spectra, analogies with finite Heisenberg 
clusters, and aspects of entanglement. It is hoped that 
the present work will motivate further experimental stud- 
ies on lateral DQDs with more than two electrons. 

Specifically it was found f Section I V Ap that, as a func- 
tion of the magnetic field, the energy spectra exhibit a 
low-energy band consisting of a group of six states, and 
that this number six is not accidental, but a consequence 
of the conservation of the total spin and of the ensuing 
spin degeneracies and supermultiplicities expressed in the 
branching diagram (described in Section IIV[) . These six 
states appear to cross at a single value of the magnetic 



field, and the crossing point gets sharper for larger inter- 
dot distances. As the strength of the Coulomb repulsion 
increases, the six states tend to become degenerate and a 
well defined energy gap separates them from the higher- 
in-energy excited states. 

The formation of the low-energy band is a consequence 
of the localization of the four electrons within each dot 
(with two electrons on each dot). The result is forma- 
tion (with increasing strength of the Coulomb repulsion) 
of a Wigner supermolecule, with the four localized elec- 
trons at the corners of a rectangular parallelogram. Using 
the spin-resolved pair-correlation functions, it was shown 
that one can map the EXD many-body wave functions 
to the spin functions associated with four localized spins 
(Section |VC]). 

This mapping led us naturally to studying analo- 
gies with finite systems described by model Heisenberg 
Hamiltonians (referred to often as finite Heisenberg clus- 
ters). Specifically, we provided a detailed interpretation 
of the EXD spin functions and EXD spectra associated 
with the low-energy band via a 4-site finite Heisenberg 
cluster characterized by two (intradot and interdot) ex- 
change integrals. More importantly, our EXD calcula- 
tions suggest a prominent oscillatory magnetic-field de- 
pendence of the two exchange integrals entering in this 
4-site Heisenberg Hamiltonian (Section IVI[) . 

Such strong independence of the exchange integrals 
has been found in previous theoretical and experimental 
studies in the simpler case of two electrons in quantum 
dots, and it has developed into a central theme in exper- 
imental efforts aiming at solid-state implementation of 
quantum computing. Our EXD results in this paper ex- 
tend the ^-dependence of the exchange integrals to larger 
numbers (N > 2) of electrons in quantum dot molecules 
(see discussion in Section [VII Al) . 

Finally, it was discussed that the EXD spin functions 
correspond to strongly entangled states known in the lit- 
erature of quantum information as A-qubit Dicke states 
(Section Eni|. 
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APPENDIX: SINGLE-PARTICLE STATES OF 
THE TWO-CENTER OSCILLATOR 

In this Appendix, we discuss briefly the energy spec- 
tra associated with the single-particle states of the two- 
center oscillator Hamiltonian given by Eq . JT]) . We follow 
here the notation presented first in Ref. |42|. For further 
details, see Ref. [43[ 

The calculated two-center oscillator single-particle 
spectrum for a double quantum dot made of two tunnel- 
coupled identical QDs (with hu> y — Tilo x i — Ulj X 2 = 3 
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FIG. 9: Single-particle spectra of a double quantum at B — 
plotted versus the distance d between two (identical) cou- 
pled QDs with a TCO confinement Uuj x i = Ttk) X 2 = hu y = 3 
meV and hi = hi = [see Eq. (JTJ)]. For all d's the bar- 
rier control parameters were taken as e\ = e\ = 0.5, i.e., 
the barrier height (depicted by the dashed line) varies as 
14(d) = Vo{d)/2. Molecular orbitals correlating the united 
(Vt = 0) and separated-dots limits are denoted along with the 
corresponding (on the right) single-QD states. Wave function 
cuts at y = along the a;-axis at several distances d (see ar- 
rows) corresponding to the lowest bonding and antibonding 
eigenvalues (solid and dashed lines, respectively) are displayed 
at the top. Energies in meV and distances in nm. 



meV) plotted versus the distance, d, between the centers 
of the two dots, is given in Fig. [5J In these calculations, 
the height of the barrier between the dots varies as a 
function of d, thus simulating reduced tunnel coupling 
between them as they are separated; we take the barrier 
control parameter as e\ — e\ = 0.5. In the calculations 
in this Appendix, we used GaAs values, m* = 0.067?7j e 
and a dielectric constant n = 12.9. For the separated 
single QDs (large d) and the unified QD (d = 0) limits, 
the spectra are the same, corresponding to that of a 2D 
harmonic oscillator (being doubly degenerate for the sep- 
arated single QDs) with a level degeneracy of 1, 2, 3, ... . 
In analogy with real molecules, the single-particle states 
in the intermediate region (d > 0) may be interpreted as 
molecular orbitals (MOs) made of linear superpositions 
of the states of the two dots comprising the DQ This 
qualitative description is intuitively appealing, though it 
is more appropriate for the weaker coupling regime (large 
d); nevertheless we continue to use it for the whole range 
of tunnel-coupling strengths between the dots, including 
the strong coupling regime where reference to the states 
of the individual dots is only approximate. Thus, for 
example, as the two dots approach each other, the low- 
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FIG. 10: Single-particle spectrum of the d = 70 nm (hu x i = 
hjjj X 2 = huiy — 3 meV, Vt, = 2.43 meV, hi = hi = 0) double 
quantum dot versus B (in T). The hui c /2 (Ml = 0, lower line) 
and Shuic/2 (Ml = 1, upper line) first and second Landau 
levels are given by the dashed lines. 



est levels (n x , n y ) with n x = n y = on the two dots 
may combine symmetrically ( "bonding" ) or antisymmet- 
rically ("antibonding") to form [0,0;0] and [0,0;1] MOs, 
with the third index denoting the total number of nodes 
of the MO along the interdot axis (x) , that is, 2n x + X, 
X = or 1; for symmetric combinations (X — 0), this 
index is even and for antisymmetric ones (I = 1), it is 
odd. Between the separated-single-QDs and the unificd- 
QD limits, the degeneracies of the individual dots' states 
are lifted, and in correlating these two limits the num- 
ber of x-nodes is conserved; for example the [0,0;1] MO 
converts in the unified-QD limit into the (1,0) state of a 
single QD, the [1,0;2] MO into the (2,0) state, and the 
[0,1;1] MO into the (1,1) state (see Fig. [SJ. Note that 
MOs of different symmetries may cross, while they do 
not if they are of the same symmetry. 

In a magnetic field, the TCO model consitutes a gener- 
alization of the Darwin- Fock model^ for non-interacting 
electrons in a single circular QD. The single-particle spec- 
tra for the DQD (d = 70 nm, V b = 2.43 meV) in a mag- 
netic field (B) are shown in Fig. [TO] (here we neglect the 
Zeeman interaction which is small for our range of B 
values with g* = —0.44 for GaAs). The main features 
are: (i) the multiple crossings (and avoided crossings) as 
B increases, (ii) the decrease of the energy gap between 
levels, occurring in pairs (such as the lowest bonding- 
antibonding pair), portraying an effective reduced tun- 
nel coupling between the QDs comprising the DQD as B 
increases, (iii) the "condensation" of the spectrum into 
the sequence of Landau levels (Ml + l/2)7kj c , Ml = 0, 
1, 2, ... (the Ml — and Ml = 1 bands are depicted, 
respectively, by the lower and upper dashed lines in Fig. 
[TO)) . This is similar to the behavior of the single-particle 
Darwin-Fock spectrum for harmonically confined elec- 
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trons in a circular QD^ (note however that the geometry harmonic confinement), 
of the DQD is non-circular and deviates from a simple 
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